Molecular Epidemiology and Antifungal Susceptibility Profile in Nakaseomyces glabrata Species Complex: A 5‐Year Countrywide Study

ABSTRACT Background The current study aimed to identify Iranian Nakaseomyces (Candida) glabrata complex species in the clinical isolates and determine their antifungal susceptibility profile. Methods In total, 320 N. glabrata clinical isolates were collected from patients hospitalized in different geographical regions of Iran. The initial screening was performed by morphological characteristics on CHROMagar Candida. Each isolate was identified by targeting the D1/D2 rDNA using a multiplex‐PCR method. To validate the mPCR method and determine genetic diversity, the ITS‐rDNA region was randomly sequenced in 40 isolates. Additionally, antifungal susceptibility was evaluated against nine antifungal agents following the CLSI M27‐A4 guidelines. Results All clinical isolates from Iran were identified as N. glabrata. The analysis of ITS‐rDNA sequence data revealed the presence of eight distinct ITS clades and 10 haplotypes among the 40 isolates of N. glabrata. The predominant clades identified were Clades VII, V, and IV, which respectively accounted for 22.5%, 17.5%, and 17.5% isolates. The widest MIC ranges were observed for voriconazole (0.016–8 μg/mL) and isavuconazole (0.016–2 μg/mL), whereas the narrowest ranges were seen with itraconazole and amphotericin B (0.25–2 μg/mL). Conclusion Haplotype diversity can be a valuable approach for studying the genetic diversity, transmission patterns, and epidemiology of the N. glabrata complex.


| Introduction
Over 150 yeast species have the potential to cause disease in humans [1].Nearly 90% of all Candida infections are caused by C. albicans, Nakaseomyces (Candida) glabrata, C. parapsilosis, and C. tropicalis [2].In recent years, non-albicans Candida species, including N. glabrata, have been considered important opportunistic pathogens and second rank in invasive candidiasis, which are emerging threats to resistance to fluconazole [3][4][5][6].N. glabrata is known as one of the leading causes of invasive fungal infections because of its high mortality rate [7,8].In addition, C. auris is an emerging fungus that presents a serious global health threat [9].
A decade ago, advancements in genomic research revealed that N. glabrata belongs to the Nakaseomyces clade of Saccharomycotina, which also includes five other species such as C. nivariensis and C. bracarensis [10][11][12].Because of several different traits between these cryptic species in their antifungal profile and mortality rate, differentiating them is very important [11].Recently, a multiplex PCR (mPCR) protocol has been used to distinguish N. glabrata from C. bracarensis and C. nivariensis, which are closely related species [13].
Although echinocandins are now considered a first-line empirical treatment for invasive Candida infections, there has been an increase in resistance to echinocandins in Candida species [14][15][16][17].In addition, considering that more than a decade has passed since the description of these cryptic species, there is still limited information about antifungal susceptibility patterns and distribution in different geographical areas.From an evolutionary standpoint, detecting and identifying these cryptic species can be useful to clarify their pathogenicity [18,19].Because of the increasing prevalence of azole resistance, the treatment of invasive N. glabrata infection is considered an important clinical challenge [20].To determine the source of fungal infection, drug susceptibility, geographical distribution, and any relationships between them, it is necessary and important to determine the genotype [21].Therefore, the present multicenter study aimed to identify clinical isolates of Iranian N. glabrata complex species, perform genetic diversity analysis, and characterize their antifungal susceptibility profiles.

| DNA Extraction
Genomic DNA from N. glabrata isolates was extracted using the phenol-chloroform DNA method [22].For this purpose, the colonies were lysed in 300 μL of lysis buffer (containing 200 mmol/L Tris-HCl [pH 7.5], 25 mmol/L EDTA, 0.5% [w/v] SDS, and 250 mmol/L NaCl), followed by incubation at 100°C for 15 min and then centrifugation.The supernatant was mixed with 20 μL of 3 M sodium acetate, frozen at −20°C for an hour, and centrifuged at 12,000 g for 10 min.The supernatants were precipitated with an equal volume of ice-cold isopropanol, centrifuged at 10,000 g for 10 min, washed with ice-cold 70% ethanol, air-dried, and resuspended in 50 μL TE buffer.A nanodrop spectrophotometer analysis was performed to assess the quality and quantity of the isolated DNA.Finally, the DNA samples were stored at −20°C until needed.

| Multiplex Polymerase Chain Reaction
Multiplex PCR was utilized to identify N. glabrata species by targeting the D1/D2 domain of the large subunit ribosomal DNA with specific primers, enabling discrimination between N. glabrata, C. bracarensis, and C. nivariensis, as previously described [23].The PCR reaction was carried out in a final volume of 25 μL as follows: 8.5 μL of master mix, 12.5 μL of ultrapure water, 1 μL each of N. glabrata complex primers, and 1 μL of DNA template.The PCR program included a predenaturation step at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 60°C for 30 s, extension at 72°C for 30 s, and a final extension at 72°C for 8 min.PCR products were separated on 2% agarose gels stained with SafeStain (Thermo Fisher Scientific, MA, USA) and subsequently visualized using gel documentation to differentiate between species complexes on the basis of the size of the amplicon.The mPCR is designed to produce amplicons of approximately 386, 214, and 588 bp from N. glabrata sensu stricto, C. bracarensis, and C. nivariensis, respectively.

| PCR Sequencing of ITS Region of rDNA
DNA sequencing was performed for the ITS-rDNA region of 40 randomly selected isolates to confirm the N. glabrata species complex and investigate the genetic diversity among the isolated species.For this purpose, the ITS region was amplified using the pan-fungal ITS1 and ITS4 primers, as described previously [24].PCR products were sequenced using the Sanger method in the forward direction.Nucleotide BLAST analysis was performed on the resulting sequences from each isolate by comparing them to sequences in the NCBI database.Accurate species identification requires a minimum sequence identity of at least 99% for reference sequences.The obtained sequences were submitted to the GenBank database.

| Phylogenetic Analysis
The clade and variation among the N. glabrata isolates were studied by comparing ITS-rDNA region sequences.The full length of sequences was cut and trimmed as required using MEGA 11 software.Multiple sequence alignments were carried out to compare the 40 randomly selected N. glabrata isolate sequences.

| ITS-rDNA Haplotype Diversity Analysis
To assess the predominant ITS haplotypes, determine intraspecies differences, and assess genetic diversity, an analysis of haplotype diversity was conducted.Pairwise comparisons and multiple sequence alignments were performed using CLUSTAL W2.For this purpose, a set of 40 randomly selected N. glabrata species were used, and their nucleotide sequences of ITS-rDNA were aligned using the MEGA 11 software.The aligned sequences were exported as a Nexus format file.Then, the NEXUS format was used in the PopART v. 1.7 software to prepare visualizations of the haplotype networking and genetic distance among Iranian N. glabrata species complex isolates [26].
To ensure quality control and reproducibility of the testing methods, American Type Culture Collection (ATCC) control strains (C.krusei ATCC 6258 and C. parapsilosis ATCC 22019) were included.The minimum inhibitory concentrations (MICs) were interpreted on the basis of the proposed CLSI M27-A4 and CLSI M60 breakpoints [27,28].

| Statistical Analysis
The data, including the ITS clade, source of N. glabrata species complex isolation, geographic location of the isolates' origin, and MIC of the isolates, were transferred to IBM SPSS Statistics software version 24.The probability of isolating a distinct ITS type with respect to the type of N. glabrata, location of isolation, and MIC were computed using the chi-squared test.An alpha level of 0.05 was used as a statistical significance criterion for all statistical tests.

| Phenotypic Identification of Nakaseomyces glabrata Complex Isolates
The purity of 320 clinical isolates was verified by their growth on CHROMagar Candida medium.Table 1 shows different clinical samples among the 320 identified N. glabrata complexes.The blood sample had the highest sample count at 23.1%, whereas CSF had the lowest count at 0.93%.

| Establishment of mPCR Assay and Analysis of Clinical Nakaseomyces glabrata Isolates
All 320 isolates produced ~386 bp amplicons, which is characteristic of N. glabrata sensu stricto strains.None of the clinical isolates produced the 588 and 214 bp amplicons that are the characteristics of C. nivariensis and C. bracarensis, respectively.The accession numbers for the 40 ITS sequences were released in the NCBI GenBank nucleotide database (https:// www.ncbi.nlm.nih.gov/ nuccore) under accession numbers OR647599-OR647638.

| Phylogenetic Analysis
On the basis of ITS region phylogenetic analysis data, eight distinct clades (arbitrarily assigned as I-VIII) were identified among 40 clinical N. glabrata isolates.Figure 1 illustrates a maximum likelihood phylogenetic tree constructed using DNA sequence data for the ITS-rDNA region from a total of 40 isolates of N. glabrata sensu stricto.On the basis of the findings from the phylogenetic tree (Figure 1), a total of eight distinct clades were identified.Among these clades, Clade VII accounted for 22.5% of the isolates, followed by Clade VI and Clade IV, each representing 17.5% of the isolates.Clades I and II accounted for 12.5% each, whereas Clade III represented 10% of the isolates.Clade VIII was observed in 5% of the isolates, and Clade V was found in 2.5% of the isolates.
Figure 2 demonstrates the genetic diversity analysis of 40 N. glabrata sensu stricto isolates, focusing on two aspects: city (A) and sample sources (B).Tehran had the highest number of clades, with a total of seven clades identified among the isolates (Figure 2A).On the contrary, the city of Ahvaz had the lowest number of clades, with only two clades observed.Clade V was exclusively isolated in Tehran and was also identified in isolates obtained from catheters.Clade V was exclusively detected in samples collected from Tehran, indicating its specific presence within this geographical location.On the contrary, Clade IV was identified in isolates from all cities under investigation, suggesting a wider distribution across different geographical regions.These findings highlight the geographical specificity of Clade V to Tehran and the broader presence of Clade IV across all cities.
Clade II was exclusively identified in 100% of BAL samples, whereas Clade VII was solely detected in 100% of urine samples (Figure 2B).Additionally, Clade III was exclusively present in 100% of CSF samples.A statistically significant relationship was detected between clade and sample type, as indicated by a p value of 0.01.Among women, Clade VII was the most prevalent, whereas Clade VI was the most common among men.No statistically significant relationship was found between clade and gender (p value = 0.4).

| ITS Haplotypes Among Nakaseomyces glabrata Isolates
Regarding the ITS sequence, 10 haplotypes were detected on the basis of the sequence analysis of 40 isolates obtained from different cities.Investigation of haplotype diversity within the strains   led to the identification of 10 haplotypes, six of which were dominant and two of which were common.The largest haplotype groups were related to Haplotypes 8 and 7, which are related to Clades I, IV, V, VI, and VIII, which belonged to the populations of Tehran, Mashhad, Ahvaz, Shiraz, and Sari (Figure 3).A network analysis was performed on the haplotypes found during the sequencing process to validate the results obtained from the genetic distance matrix method.
The haplotyping analysis demonstrated a nucleotide diversity of 0.0025 between different strains of N. glabrata.Generally, the results indicated a high rate of genetic diversity and haplotype diversity.Analysis of ITS sequences revealed that both indices (π = 0.0025) were high.Tests of the standard neutral model for a demographically stable population-Tajima's D (0.730) were carried out on the ITS sequences.The neutrality statistics were negative and not significant at p = 0.47, perhaps due to the intronic regions present in the ITS sequences.

| Antifungal Susceptibility Testing
All isolates included in the study exhibited sensitivity to echinocandins.Three isolates belonging to Clades I and VII exhibited high Minimum Inhibitory Concentration (MIC) values for fluconazole.Among all the clades, fluconazole was the only drug that demonstrated relatively high MIC values.

| Discussion
Nakaseomyces glabrata is widely recognized as one of the primary causes of invasive candidiasis, a condition that is becoming an increasingly significant public health concern globally.Two cryptic species, C. bracarensis and C. nivariensis, are not distinguishable from N. glabrata [11,29].Molecular identification and differentiation between these species' complexes are crucial because of variations in antifungal susceptibility and virulence profiles [1,18,30].Because of the differences in mortality rate and antifungal susceptibility among N. glabrata complex species, as well as the lack of comprehensive epidemiological data and conflicting findings in various studies, the antifungal susceptibility profile and geographic distribution of these cryptic species are important.Therefore, accurate identification of the species in this complex is necessary for proper antifungal treatment [1,31,32].In relation to this matter, techniques such as multiplex PCR enable quick, affordable, precise, and reliable identification of cryptic species of N. glabrata.In spite of C. nivariensis and C. bracarensis reports from different countries, limited information has been provided about their antifungal susceptibility patterns, epidemiology, actual prevalence, and biological niches.Therefore, the present study was conducted with the aim of identifying and determining the clinical characteristics, antifungal susceptibility, and clade and haplotype of the N. glabrata species complex in Iran.In the present study, 320 isolates collected from various regions of Iran were analyzed, and none of them were identified as C. bracarensis or C. nivariensis.These findings align with recent research, in which similarly, Esposto et al. [33] investigated 1000 clinical isolates obtained from 14 regions in Italy but did not detect any instances of C. nivariensis or C. bracarensis.Furthermore, in a study conducted in Kuwait, among 440 clinical isolates, 100% of the isolates were reported as N. glabrata sensu stricto [34].In a study of 143 clinical isolates collected in Spain, three samples (2%) were identified as C. bracarensis, whereas no isolates were identified as C. nivariensis [35].Additionally, Angoulvant, Guitard and Hennequin.[8] reported the prevalence of these two cryptic species as NC.bracarensis (0.01%) and C. nivariensis (0.12%).Various molecular techniques have been previously documented for the Hp-4 3 ( identification of isolates specific to N. glabrata sensu stricto, C. nivariensis, and/or C. bracarensis [13,30,36,37].Among these methods, similar to our study, Arastehfar et al. [23] and Asadzadeh et al. [34] utilized a simple, fast, and cost-effective method on the basis of mPCR for the accurate identification of the N. glabrata complex, which we confirm as well.Moreover, according to the present study, no isolates of C. nivariensis and C. bracarensis were isolated, indicating the low prevalence of these cryptic species in Iran.
Given that there is considerable variation in virulence patterns among N. glabrata strains, as well as associations between specific clades and increased mortality rates, the critical role of genotyping techniques in diagnosis has been highlighted [38,39].
Considering that few studies have been done on the genotyping of the N. glabrata complex and there is little information about it, two studies close to ours can be mentioned in this regard.Asadzadeh et al. [34] in Kuwait, 85 isolates of N. glabrata were randomly selected from among 440 isolates using ITS and sequencing.Finally, 28 haplotypes were identified, of which 18 isolates had distinct clades and 67 isolates were grouped into 10 haplotypes.In this study, 40 samples were randomly selected from a total of 320 N. glabrata isolates, resulting in the identification of eight clades and nine haplotypes.The variation in the number of clades and haplotypes observed in different studies may be attributed to factors such as the number of samples, duration of sample collection, and inclusion of samples from diverse sources.These findings support the high degree of heterogeneity among the N. glabrata species complex in Iran (Tehran, Shiraz, Yasuj, Isfahan, Ahvaz, Mashhad, and Mazandaran provinces), despite the same source of samples and geographical locations, suggesting that isolates from patients originated from another side during hospitalization.
The strength of our study is that it is the first investigation of the clade of the Iranian species complex of N. glabrata from Tehran, Shiraz, Ahvaz, Mashhad, and Mazandaran provinces.
In order to assess the genetic diversity and transmission patterns of the N. glabrata species complex, genotyping is essential for epidemiological investigations, in particular where patients are being treated with different treatments.The limitation of this study is limited to the number of isolates analyzed and a lack of information about ITS sequences for N. glabrata species collected in Iran.
There is not sufficient data to support the hypothesis that all species of the N. glabrata complex are generally susceptible to echinocandins.In light of different susceptibility patterns among the N. glabrata species complex, the lower activity of echinocandin was reported against C. nivariensis [40].On the basis of the findings of different studies and various drug susceptibility patterns of the N. glabrata complex, accurate identification of these species is very important [41].According to a study conducted by Hou et al. [18] using molecular methods, 12 isolates of C. nivariensis and 1 isolate of C. bracarensis were identified, and 100% of them were susceptible to echinocandins by antifungal susceptibility evaluation.According to the findings of this study, N. glabrata sensu stricto isolates showed resistance to caspofungin and micafungin while being sensitive to anidulafungin.These findings are consistent with previous studies.It is noteworthy that recent evidence casts doubt on the future effectiveness of echinocandins, which are a class of antifungal drugs.Several studies have reported resistance to azoles in these two species [1,30,[40][41][42].The most striking result in our study was that all N. glabrata had a low MIC value for itraconazole, voriconazole, and posaconazole; all the isolates were found to be susceptible to echinocandins, whereas 7.5% of the samples demonstrated resistance to fluconazole.According to various reports, as there is a significant difference between these two species in terms of susceptibility to antifungal drugs compared with N. glabrata senso stricto, studies with larger sample sizes in different geographic locations should be conducted to clarify their epidemiology and antifungal susceptibility profile [30].Our ability to identify and share data on a national and global scale will help us learn about these species and improve diagnosis and treatment.It can also aid in understanding the evolution and spread of this fungal pathogen.

| Conclusions
Detection and identification of cryptic N. glabrata species in relation to their closely related counterparts present challenges that have led to uncertainties concerning their susceptibility to antifungal agents and their epidemiology.In addition, because of the importance of identifying these species with the shortest turnaround time, it is very necessary to develop cost-effective, efficient, and reliable methods.Given the documented increase of these cryptic species in multiple European countries and their propensity to develop resistance to antifungal treatments, it is recommended to continue monitoring these emerging cryptic species.

FIGURE 1 |GU199443. 1 GU199441. 1 GU199440. 1 MW284463. 1
FIGURE 1 | Maximum likelihood phylogeny tree based on DNA sequence data for the ITS region of rDNA from 40 Nakaseomyces glabrata sensu stricto isolates.The associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches.

FIGURE 2 |
FIGURE 2 | Genetic diversity of 40 Nakaseomyces glabrata sensu stricto isolates based on city (A) and sample sources (B).

FIGURE 3 |
FIGURE 3 | Haplotype network generated for Internal transcribed spacer (ITS) region sequences on 40 Nakaseomyces glabrata sensu stricto isolates.The size of the circle indicates the relative frequency of sequences belonging to a particular haplotype (smallest circle = 1 sequence to largest circle = 11 sequences).Hatch markers along the network branches indicate the number of SNP.

TABLE 2 |
In vitro susceptibilities of nine antifungal drugs against 40 clinical strains of Nakaseomyces glabrata.